Optimization Schemes for Selective Molecular Cleavage 
with Tailored Ultrashort Laser Pulses 



Kevin Krieger*^, Alberto Castro'^'*, E. K. U. Gross^ 

"■Max-Planck-Institut fiir Mikrostrukturphysik, Weinberg 2, D-06120 Halle, Germany. 
^Institute for Biocomputation and Physics of Complex Systems (BIFI), University of 
Zaragoza, E-50018 Zaragoza, Spain. 



Abstract 

We present some approaches to the computation of ultra-fast laser pulses 
capable of selectively breaking molecular bonds. The calculations are based 
on a mixed quantum-classical description: The electrons are treated quan- 
tum mechanically (making use of time- dependent density- functional theory), 
whereas the nuclei are treated classically. The temporal shape of the pulses 
is tailored to maximise a control target functional which is designed to pro- 
duce the desired molecular cleavage. The precise definition of this functional 
is a crucial ingredient: we explore expressions based on the forces, on the 
momenta and on the velocities of the nuclei. The algorithm used to find the 
optimum pulse is also relevant; we test both direct gradient-free algorithms, 
as well as schemes based on formal optimal control theory. The tests are 
performed both on one dimensional models of atomic chains, and on first- 
principles descriptions of molecules. 
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1. Introduction 

Soon after its first operation, the laser was expected to become the 
ultimate surgical tool at the nanoscopic level: Light, at convenient wave- 
lengths, monochromatic, coherent, and intense, [2;] was believed to open the 
avenue to selectively break (or create) molecular bonds. Unfortunately, the 
early attempts to perform this kind of photo-chemistry were only occasionally 
successful. 3|, 14] These attempts used "simple" monochromatic lasers, tinker- 
ing only with two parameters: the frequency and the intensity. However, the 
energy, tuned to a particular vibrational frequency and initially deposited on 
the corresponding bond, is soon re-distributed to the rest of the modes, and 
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produces undesired global heating instead of selective cleavage. [5[ 

The "controlled" laser assisted photo-chemistry advanced along with im- 
provements on laser technology, with methodologies such as the control of 



quantum interference proposed by Brumer and Shap iro, 



,_,17|,18J the "pump- 
dump" control proposed by Tannor and Rice. 9 . id l stimulated Raman adi- 
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abatic passage, [Uj wave-packet interferometry. [12l| and others. [13l. |lj] The 
key ingredients, beyond mere mono-chromaticity and intensity, were shown 
to be coherence (and therefore, interference), detailed shaping, and ultra- 
short pulse duration (in the femto-second time scale). The most successful 
technique is adaptive feedback control (AFC), as proposed by Judson and 
Rabitz. [isl and first realised in 1997. 16] 



There are two important components in an AFC experiment: the pulse 
shaper,|l3] and the search algorithm fed by the repeated measurement out- 
come. The former is an instrument that allows to almost arbitrarily design 



laser pulses. The increasing versatility of modern laser sources (regarding 
pulse length, power, and accessible frequencies), and the capacity of pulse 
shapers to modify the produced pulses, set the boundaries that theoretical 
studies such as the one presented in this work must respect; however these 
boundaries are rapidly pushed further, allowing more versatile pulses. 



Quantum optimal control theory (QOCT) |l8l . Il9l |20| . l2lj is the most gen- 
eral theoretical framework aimed to the prediction of laser pulses that are 
optimal for a given task. It is the translation to the quantum realm of a 
very broad mathematical area, optimal control, that is best formulated in 



the language of systems theory, 
initiated in the 80s 
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23| Its use for quantum processes was 



2l| - responding to the initial experimental stir. 



In some way, QOCT encompasses all the previously mentioned optimisation 
methods (inasmuch as it may describe them theoretically). The theory is 
constructed on top of some chosen level of approximation for the descrip- 
tion of the process that is to be optimised. Here lies the main limitation of 
QOCTijl^ it may only be predictive if the system is simple enough to allow 
for an accurate approximation of its evolution. In most cases, however, the 
process is too complex. 

If some reliable predictive power is to be expected from any QOCT calcu- 
lation, one should attempt a first-principles description. In particular, in the 
regime of interest, the dynamics of the electrons should be carefully treated: 
high intensity electric fields at high frequencies affect directly the electronic 
degrees of freedom. Indeed, when many-electron systems are irradiated with 
strong femtosecond pulses a number of interesting non-trivial photo-reactions 
may take place: above-threshold or tunnel ionisation, bond hardening or soft- 
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ening, high harmonic generation, photo-isomerisation, photo-fragmentation, 



Coulomb explosion, etc. 
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26| . Yet most of the computational work un- 



til now has relied on simplified models, and has usually worked with nuclear 
wave packets - defined on a few relevant reaction variables, after a reduc- 
tion of dimensionality has been postulated - moving on one or a few Born- 
Oppenheimer potential energy surfaces, and therefore mostly ignoring the 
dynamic behaviour of the electrons. Direct, first-principles, electronic con 

one-electron cases. 



28 



29 



3Q| 



One viable alternative to treat electronic motion in an ab initio way is 
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32| Recently some 



33| 



time-dependent density-functional theory (TDDFT) 
of us have demonstrated the feasibility of performing QOCT with TDDFT 
This was not obvious due to the non-linear character of the TDDFT equa- 
tions: the usual QOCT equations assume a standard, linear Schrodinger-like 
evolution, and the resulting QOCT equations are correspondingly simple. 
However, the presence of the Hartree, exchange and correlation term in the 
TDDFT equations need special care. 

TDDFT offers reasonable accuracy when dealing with the non-linear re- 
sponse of molecular systems, with a fraction of the cost of methods based 
on the wave function. Furthermore, the electronic system described within 



TDDFT may then 
classical description. 



3e coupled to the ionic motion in a mixed quantum- 
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This model will obviously ignore quantum 
nuclear effects, but may be sufficient for the description of many processes. 
In this work we present our first results based on this combination. In Sec- 
tion m the essential equations are displayed, as well as a brief description of 
the numerical procedure. Section [3] describes the results of the optimisations 
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when the target functional is defined in terms of the values of the forces on 
the nuclei at the end of the laser pulse, for ID models, whereas in Section HI 
the target functional is defined in terms of the momenta. In Section [5|, the 
attempt to selectively break molecular chains is described. Finally, Sections 
M and [7] display results for fully ab initio 3D calculations. 

2. Methodology 

2.1. Essentials of QOCT 

We consider a quantum mechanical system governed by Schrodinger's 
equation during the time interval [0,T] (atomic units will be used hereafter): 

= H[u,t]^{x,t) , (1) 
*(x,0) = ^oix), (2) 

where x is the full set of quantum coordinates, and -u is a control, typically a 
set of parameters that determine the precise shape of an external potential 
applied to the system. Mathematically, we can distinguish two types of 
"representation" for the control u: 

1. u is a real valued continuous function defined on the time interval of 
interest (the control function); we will call this a "real-time" repre- 
sentation of the control. For example, the Hamiltonian may have the 
form: 

H[u,t] = Ho + u{t)D . (3) 

2. M is a set of real parameters that modifies the precise shape of the 
Hamiltonian; typically, this set of parameters fixes the form of a con- 
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trol function; we will call this a "parameterised" representation of the 
control. 

In any case, the specification of w, together with an initial value condition, 
^(0) = determines the full evolution of the system, ^[u], via the propa- 
gation of Schrodinger's equation. 

We wish to maximize the function G, 

G[u]^FMu],u], (4) 

where F is the so-called "target functional" ; in many cases it is spht into two 
parts, = + J2[u\, so that Ji only depends on the state of the 

system, and J2 is called the "penalty" , and depends explicitly on the control 
u. An important distinction should be made regarding Ji. 

1. It may depend on the full evolution of the system during the time 
interval [0,T]; this is usually called a time-dependent target. We may 
write this as Ji[^] = Jj' [^'], where the J{ ' [^] functional admits 

continuous functional derivatives, in particular ] is continuous 

d^/*{x, t) 

at t = T. 

2. Ji may only depend on the state of the system at the end of the prop- 
agation, which we may write as = Ji[if{T)]. 

Of course, Ji may be defined as a combination of the two options, i.e.: 

= jFi[*]+jn*(T)]. (5) 

Note that, in this case: 

" (6) 
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In most cases these functionals are defined as the expectation value of 
some observable O. For example: 

jn^(r)] = (M/(T)|6|M/(T)), or: 



(7) 



Jo 

One needs now an optimization algorithm to find the maximum (or max- 
ima) of G. Two broad families can be distinguished: gradient-free procedures, 
that only require some means to compute the value of G given a control input 
u, and gradient-based procedures, that also necessitate the computation of 
the gradient of G with respect to u (more precisely, the functional derivative 



T 



if -u is a continuous function in time). We wi 
that can be found elsewhere in several forms; 
equations are: 



1 not repeat here a derivation 
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39j the key 



2Im / dt {x[u]it)\VuH[u,t]mum) 



(9) 



in case -u is a set of real parameters, and: 

6G 



+ 



Su{t) 5u{t) 

2\m{x[u]{t)\D\m[u]{t)) , (10) 

if is a function in time, and the Hamiltonian is given by Eq. |3l 

Note that a new "wave function", xVAi has been introduced; it is given 
by the solution of: 



x[u]{x,T) 



T 
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(11) 
(12) 



This is similar to the original Schrodinger's equation (Eqs. [Hand [2]), except: 



(1) It may 
dependent 



\C\ T'] 

De inhomogeneous, if J{ ' is not zero (i.e. if the target is time- 



38, 



39|), and (2) The initial condition is given at the final time 
t = T, which implies it must be propagated backwards. 

The computation of the gradient or functional derivative of G, there- 
fore, requires \E'['u] and which are obtained by first propagating Eq. [1] 
forwards, and then Eq. [11] backwards. The maxima of G are found at the 
critical points VnG'['u] = or = 0; in order to arrive to these maxima 
one can use a variety of algorithms, some of which are listed in Section 12.3.21 

2.2. Mixed quantum- classical description with TDDFT 

Instead of solving the many-electron Schrodinger equation, TDDFT al- 
lows to work with a set of one-electron equations, the Kohn-Sham (KS) 
system, corresponding to a fictitious system whose one-particle density is by 
construction identical to that of the real one: 

= -^VVi(r,i) + Kxt(r,t) 

+ ^^HartreeK](^ + V^c[n]{f, t)] ipi{f, t) , (13) 
N 

n{f,t) = ^2|^,(f,t)|2 = nt(f). (14) 

1=1 

We will assume a system with 2N electrons in a spin compensated configu- 
ration, evolving in a spin independent Hamiltonian. This means doubly 
occupied KS orbitals Lpi, i = 1,. . . ,N. The system evolves on an external 
time-dependent potential Vext, that may include the interaction with a set 
of nuclei, as well as external electric fields. The Hartree term fHartree is the 
classic electrostatic potential, and the rest of the electron-electron interac- 
tion is encoded in the exchange and correlation potential Wxc- In this work, 
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we will only use the adiabatic extension of the local density approximation 
(LDA) , 40| although the extension to other more sophisticated schemes is 
straightforward. 

The external potential can depend on a control function, and therefore 
control theory can be employed to find optimal evolutions of the KS sys- 
tem. Note, however, that the KS equations are not akin to the conventional 
Schrodinger equation, since they are non-linear. The QOCT expressions 
derived above are therefore not valid; the correct equations have been pre- 
sented elsewhere; 



331 ] however, in this work we will either (1) take the inde- 
pendent electron approximation, which amounts to ignoring the mentioned 
non-linearity, for model calculations, or (2) utilize a gradient-free version of 
QOCT, for which we can use the full-fledged version of TDDFT. 

In order to describe the combined coupled movement of electrons anc 
(classical) nuclei, one can perform Ehrenfest dynamics on top of TDDFT. jsT] 
The external term Vext will couple the electrons to N^^c nuclei located at 
positions Ra{t) through an expression in the form: 

VeAr,t; {Mt)}) = J2 , g ^ + m ■ r (15) 

The evolution of the nuclear positions is then governed by an Ehrenfest equa- 
tion in the form: 

U \Rc{t)-RM' 

- y"dV n(f, t)V^^v,xtir, t; {Rpit)}) . (16) 



9 



2.3. Numerical implementation 

All the ideas described above have been implemented in the octopus 
code. Since the numerical details of this platform are described elsewhere 
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42| here we will only list some essential points. The laser field and the 



optimization algorithms are described below with more detail. 

• Wave functions and densities are represented on a regular rectangular 
real space mesh. This is a suitable scheme to describe high intensity 
laser-electron interactions, since the electronic density visits regions in 
space far from the localized basis sets typically used in other schemes. 
Furthermore, the intrinsic locality allows for easy parallelisation, and 
the only parameters controlling convergence are the grid spacing and 
the simulation box size. 

• The electron-ion interaction is modelled with pseudopotentials. In this 
way, the Coulomb singularity is avoided, and the core electrons are 
removed from the calculation. For the first results described below, 
however, we will use ID models, and the soft-Coulomb interaction to 
avoid singularities. 

ved in real time with the help of a number of 



43| This is crucial since all algorithms require 



The KS orbitals are evo 
propagating algorithms, 
multiple propagations. 

The code performs realistic 3D calculations, but it also allows ID and 
2D models, such as the ones we will present below. 
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2.3.1. The laser field. 

We will assume that laser pulses can be described in the dipole approx- 
imation, which is valid given the wave lengths and intensities that will be 
considered. In consequence, it suffices with an electric field in the form: 

m = eit)p, (17) 

where p is a unit vector that determines the polarization direction, and e{t) 
determines the temporal dependence, and is the object to be optimized - i.e. 
the control function. 

Not any function in time is admissible as a solution; there are physical and 
experimental constraints that must be respected. For example, an important 
physical constraint is: 

T 

dt e{t) = 0. (18) 

This condition follows from Maxwell's equations for a freely propagating 
pulse in the electric dipole approximation. [44'| Also, the pulses must obviously 
start and end at zero: 

e(0) = 6(T) = . (19) 

It is important to reduce the search space to functions that are experimentally 
accessible, which means a limitation on the accessible frequency components, 
and on the intensities. Regarding the latter, usually it is done by considering 
the integrated intensity or fluence, defined as: 

J^[e]= fdte^t). (20) 



Spectral constraints can also be imposed either by penalizing the undesired 
frequencies in the definition of the target, or by restricting from the start 
the search space to the correct subspace. 
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As discussed earlier, we may use a real-time representation, and therefore 
e{t) is directly the control object u, or a parameterised representation, in 
which this control function e{t) is determined by a set of parameters u. This 
distinction is relevant for the mathematical derivations (since in the former 
case functional derivatives must be used, whereas in the latter case one uses 
normal gradients). Numerically, however, a function in real time must also be 
discretized, and therefore the distinction disappears. Nevertheless, typically 
the number of degrees of freedom (number of grid points in time) will be 
much larger, and therefore the algorithms utilized will differ. 

Regarding the choices for the parameterisation, it is a natural choice to 
expand the control field in a basis set, and to establish the coefficients of this 
expansion as the parameters: 



N is the dimension of the real basis set {gn{t)}- It is chosen to be orthonormal 
over the interval [0,T]: 



N 




(21) 



n=l 




(22) 



In our calculations, two basis sets have been used: a sine basis: 





(23) 



or a normal Fourier basis: 



9n{t) 




|cos(^nt) 



|sin(^(n 



f)t), n^{f + l),...,N. 



n 



= 1 ^ 



(24) 
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The representation in these basis sets has the advantage that spectral con- 
straints can be automatically enforced: the maximum frequency is given by 
the choice of A^, and we we will not include the zero- frequency component, 
in order to satisfy condition flTSj) . 

We can directly choose the basis set expansion coefficients as constrol 
parameters, or else constrain further the search space to meet other physical 
or experimental requirements, by defining the coefficients as functions of a 
reduced set of parameters: e„ = CnM- Our choices have been the following: 

• A constrained sine series. The sine series, Eq. [2^ automatically fulfills 
the condition given by Eq. f lT9|) . To meet condition f|T8|) . however, the 
following relation would have to be fulfilled: 

N/2-1 _ 

y -^^^ = 0. (25) 

For some of the cases presented below, we also enforced a fixed fiuence. 
As function of any orthonormal basis set coefficients, the fiuence is 
given by: 

TV 

m = Y.'n- (26) 

n=l 

Setting the fiuence to a predefined value J-q amounts to requiring the 
vector e to belong to a hypersphere. We may then transform e into 
hyperspherical coordinates; the A^ — 1 angles 6j will span the new search 
space, of one dimension less. 

• A constrained Fourier series. If the zero-th frequency is left out, a 
Fourier series, Eq. [2H automatically fulfills condition f lTSj) . Condition 
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(pra} is met if: 



N/2 

h = - J]?n- 

n=2 



(27) 



To fulfill this condition, a first parameter transformation can be defined 

by 



N/2-1 

h =: - ctn , 



e(n+i) 



n=l 



n = l,...,(iV-l). 



(28) 



In terms of the new coordinates, it is trivial to see that the fiuence is 
given by a bilinear expression: 



J^[a\ = oP" Sex . 



(29) 



for a (A^ — 1) X (A^ — 1) symmetric matrix S. It can be diagonalized 
by performing a new change of coordinates based on an orthonormal 
matrix U: 



( 



U^SU 



Si 







V 



(30) 



s(Af-i) y 

and if we now define a final change of coordinates in the form: 



/3 = LU^a , 



(31) 



where: 



L :-- 



(32) 
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then the fluence has the simple form: 

N-l 

m = Y.Pl (33) 

n=l 

Once we have this form, in order to fix the fluence to a predefined 
value J-Q one can once again make a coordinate transformation to hy- 
perspherical coordinates, and use the N — 2 angles as search space. 

2.3.2. Optimization algorithms. 

There are two broad families of optimization algorithms: gradient-free 
and gradient-based schemes. We will utilize both in the application presented 
below. 

Gradient- free. . In experimental control experiments, the gradient of the merit 
function is seldom available, and the most used gradient-free algorithms be- 
long to the "evolutionary" or "genetic" families. These are specifically de 
signed for search spaces with large number of dimensions, typically discrete. 



47| 



461, 



However, in our code we have opted for two different schemes, which 
are sufficient for a moderate number of continuous degrees of freedom: the 



48] and Powell's NEWUOA 



classic sim ple x algorithm of Nelder and Mead, 
algorithm. |49| newer and more efficient. 

Gradient-has ed. . If the control function is described in any parameterised 
representation, then we have used a standard conjugate gradient algorithm, 



50| 



the Broyden-Fletcher-Goldfarb-Shanno variant. 

However, if the control function is represented directly in real time (which 
usually implies a large number of degrees of freedom), a number of different 
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algorithms that were specifically developed within the field of QOCT (or 
adapted to it) have been proposed. These can provide very fast convergence, 
if they are applicable. In particular, very succesful techniques are the Krotov 
method 5l| and the monotonically convergente techniques proposed by Zhu 



and collaborators. 
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53| In some of the examples given below, we will use 



one of these latter techniques. [52] 



3. Results: Control targets based on the forces 

During the breaking of a bond, the forces that act on the two separating 
nuclei should have more or less opposite directions in space, i.e., a naive but 
reasonable attempt to define a bond-breaking target is to do it in terms of the 
forces: one can attempt the maximization of the force difference between the 
nuclei that must be separated, and the minimization of the forces between 
the nuclei remaining in each fragment. In this section we describe a first 
attempt in which the target includes the value of the forces only at the end 
of the action of the pulse - it is, therefore, a static target. 

We will assume that the laser pulse is short, so that during its action the 
nuclei do not move significantly; therefore the optimization calculations will 
be performed with frozen nuclei. The idea is that the pulse should be able 
to place the electrons on a dissociating state. Later, the optimized pulse will 
be tested without the fixed nuclei restriction. In this "bond breaking test 
run", therefore, the calculation was based on the mixed quantum-classical 
description described earlier. 

We used a simple ID model of a triatomic Hydrogen molecule (see Fig.[T]), 
with two non-interacting electrons. The electron-nucleus interaction is mod- 
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Figure 1: Sketch of the ID test model. The direction of the arrows indicates the direction 
of the force optimization. 

elled with a soft Coulomb potential: 

= - V(.x-l.p + r 

where x is the electronic coordinate, and Xi is the nuclear position of nucleus 
i. Since we have two independent electrons evolving in a spin- independent 
Hamiltonian, we can assume the system to be permanently in a singlet state: 
both the two electrons occupy the same orbital ^, which is initially the 
ground state. It evolves governed by the Hamiltonian: 

1 ^ 
H[e,t] = --dl + xe{t) + ^Vr,nc{x,Xi). (35) 

i=l 

The temporal dependece of the laser field is determined by the function e(t) , 
for which wc will consider in this real time representation. 

The target functional F will be divided into the object that truly needs 
to be optimized, Ji, and a penalty function J2: 

F[*,e] = Ji[*] + J2[e]. (36) 
17 



The task of J2 is to prevent unphysically large fluences: 

J^[e] = -aj^[e] = -a I dt e\t) . (37) 

Jo 

The constant a is the "penalty factor"; it is positive, and it regulates the 
weight that is put in the low fiuence condition. 
The definition that we choose for Ji is: 

M^] = {F2[^{T)] - F3[^(T)]) - \F,[^{T)] - F,[^{T)]\' , (38) 

where Fj[^(T)] is the force acting on nucleus i at the end of the pulse action, 
and is given by: 



/ ^^^^^ rY> rv* O 

This definition of J\ attempts to maximize the force difference between nu- 
cleus 2 and 3, and minimize the force between nucleus 1 and 2. There are 
some parameters in this expressions that one can experiment with: the sec- 
ond term in the right hand side of Eq. [38] could be multiplied by a weighting 
factor, or the square could be eliminated or changed by other exponent. 
Note that this type of force target is an explicit functional of the density 
n{x,T) = 2|\E'(x,T)p- this is not so relevant in the independent electrons 
approximation taken in this case, but it is in the Kohn-Sham case that will 
be discussed later. 

We must now adapt the QOCT equations (fB. (|2|). ([TT]). (1T2|) and (fTOl) to 
this particular case. Schrodinger's equation, together with its initial condi- 
tion, dl]) and ([2]), obviously do not change. The evolution equation for the 
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auxiliary x wave function is in this case given by: 



ix,t) = H[e,t]x[e]{x,t) , (40) 
x[e]ix,T) = Oix)^ix,T), (41) 



where 



0{x) = d:r:[Vnuc{x,X2) - Vnnc{x,X3)] 

-2[Fi[^(T)] - F2[^{T)]]d,[v^U^,xi)-Vr,Ux,X2)]. (42) 
Finally, Eq. ( ITOl) takes now the form: 

^ = -2ae{t) + 2lm{x[e]{t)\xmm) ■ (43) 

At the maxima, this functional derivative is null, and therefore the solution 
field will be given by: 

e(t) = -Im(x[e](i)|x|^[6](t)). (44) 
a 

In order to solve these equations, we chose the algorithm of Zhu and 



Rabitz. 531] This is a strictly monotonically convergent algorithm, as long as 
the target functional has the form of an expectation value, a condition that 
does not hold in our case. The algorithm requires an initial guess, which is 
then iteratively improved; we chose a sine wave with sine-shaped envelope 
(see Fig. [2]): 

e(°)(i) = ^osin(^t)sin(a;ot)- (45) 

We used a laser pulse duration of T = 400a. u. and an amphtude of = 7 • 
lO^^a.u.. We tested several frequencies for the initial field: = (1, 2, . . . , 9) • 
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Figure 2: Left: Convergence plot of the optimization run. Right: Initial and optimized 
laser pulses. The optimized pulse corresponds to iteration step 30. 
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Figure 3; Left: Forces on each nucleus during the bond breaking test run. The vertical 
black line indicates the end of the laser pulse. Right: Position of the nuclei during the 
bond breaking test run. 
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10~^a.u. (note that the final yield will depend on the choice of the initial 
guess) . 

We found that convergence is by no means guaranteed - only 3 of the 9 
optimization runs showed a convergent behaviour. Furthermore we observe 
that the convergence is not monotonic. This is demonstrated in Fig. [2] where 
we show, on the left panel, the convergence history for the case uq = 4 ■ 
10~^ a.u. (all other cases were qualitatively similar). The right panel shows 
the initial and the converged laser pulse. 

We use the latter to check whether or not the bond breaks; we let evolve 
the system for 1000 a.u. (i.e., also after the pulse vanishes) with moving 
nuclei. Fig. [3] displays the forces and positions of the three nuclei during this 
process. We first observe that the forces obtained in the optimization run are 
not identical to the forces computed during this bond-breaking test run, since 
in this case the nuclei have been free to move during the short laser pulse. 
However, the differences were small, which validated (for this particular case) 
our static nuclei approximation. A second important observation is that the 
amplitudes of the force oscillations before and after the end of the pulse were 
of the order of, and even larger than, the optimized forces at the end of the 
pulse. 

In Fig. [3] (right), we observe that we got a complete atomization of the 
test model in this run - which is not the objective. This negative result was 
typical of all runs: Either the test model was still bound and the nuclei just 
oscillated around their equilibrium positions for t > 400 a.u., or we got full 
atomization, as in the case presented. This latter case was triggered by a 
strong electronic ionization. 
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In view of the strong force oscillations observed, we may conclude that 
the main reason for this negative outcome is the time-independent character 
of our control target: the forces have a strong oscillatory character, and con- 
trolling them at a single moment in time does not suffice. This consideration 
leads naturally to the subject of the next section: the definition of the control 
targets in terms of the full history of the forces - their integrated values, or 
in other words, the momenta. 

4. Results: Control targets bcised on the momenta. 

In this section, we explore the option of defining the target functional in 
terms of the momenta of the nuclei at the end of the pulse. For this purpose, 
we used the same ID model defined in the previous chapter. 

The momenta are nothing else than the integrated forces: 



Jo 

and the definition of the target functional F is simply done by replacing 
forces by momenta: 



Qualitatively, however, the problem changes, since Pi[^] are functionals of 
the full evolution of the system, i.e. we confront a time-dependent target. 
The three cases presented below differ in the manner in which the laser field 
is defined or restricted, and on the optimization algorithm. 

4.1. Gradient free optimization algorithm with fixed nuclei 

In this first case, we used a parameterised representation for the control 
function (the electric field), in particular the constrained sine series (see 




(46) 



(47) 
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Section [2.3. ip : the search space is spanned by a set of hyper spherical angles 
9 = {0j}, and therefore the fluence is constant (making unnecessary the 
introduction of a penalty function J2). 

We test now a gradient-free procedure for the maximization of the func- 
tion G[9] = F[if[9],6] = Jj[^[9]], in particular the "downhill simplex" method 
from Nelder and Mead. (48] Each function evaluation amounts to one forward 
propagation (the backwards propagations are in this case unnecessary). As 
in the previous section, we used very short pulses and assumed the fixed- 
nuclei approximation during the pulse action. The optimization runs were 
followed by the corresponding "bond-breaking test runs" , in which the nuclei 
are allowed to move to check that the molecule breaks in the intended way. 

As an initial guess for the pulse, we used, once again: 



We performed several calculations with varying values of ojq: ojq = {4 . . . 19) • 
10^^ a. u.. The amplitude Aq is adjusted so that all optimizations are per- 
formed with the same (constant) fluence. The propagating time was chosen 
to be T = 200 a. u.. All functions were then expanded in a sine Fourier se- 
ries, with frequecies Un = -^n for n = 1,2, ...12. This means 11 degrees 
of freedom for the search space, once the transformation to hyperspherical 
coordinates was done. 

All optimizations converged, and of those, half of them led to the seeked 
bond destruction. We display results for one of the runs (wq = 19 ■ 10~^ a.u.), 
since all of them were qualitatively similar. Figure H] (left) shows the conver- 
gence history of F. The right side shows the initial and the optimized laser 
pulse (iteration step 100). It is clearly visible that this optimized pulse does 




(48) 
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Figure 4: Left: Convergence plot of the gradient-free optimization. Right Initial and 
optimized laser pulses. 
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Figure 5: Left: Momentum of each nucleus during the bond breaking test run. Right: 
Positions of the nuclei during the bond breaking test run. 
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not contain very high frequency components, compared to the pulse obtained 
in the forces-based optimization. This is due to the natural frequency cut-off 
imposed by the parameterisation. 

The left plot of figure [5] displays the momenta of the nuclei during the 
bond breaking test run. It is noteworthy that the momenta did not signifi- 
cantly oscillate for t < 200 observed for the forces. The right plot, in 
turn, shows the coordinates of the nuclei during the bond breaking test run. 
It is clear that the intended goal was achieved: nucleus 3 dissociates from 
nucleus 1 and 2, that stay bound. 

Despite the successes, the nuclear movement was not completely negli- 
gible during the action of the laser pulse. This can already be seen in the 
right panel of Fig. O In order to further study the influence of the nuclear 
movement, we performed runs with different pulse durations (T = 100 a. u. 
and T = 400 a.u.). The result is that for 100 a.u. many runs succeeded, while 
for 400 a.u. no run did. We may conclude that (1) constructing the control 
target functional in terms of the momenta of the nuclei is an appropriate 
approach to the problem of selective bond cleavage, but (2) the movement 
of the nuclei is, in general, not negligible when performing the optimization, 
unless the laser pulses are very short. 

4-2. Gradient free optimization algorithm with moving nuclei 

The natural next step is therefore to include the ionic motion during the 
optimization runs, in order to allow for larger pulse durations. We have 
attempted this using exactly the same model and target definition as in the 
previous section. The only difference is that, during the action of the pulse, 
the dynamic variables include not only the electronic orbital, but also the 
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nuclear coordinates and momenta. 

The laser pulse was represented in the same way as in the previous sec- 
tion: the set of hyperspherical angles that describe the fixed-norm (i.e. fixed 
fiuence) coefficients of a sine series expansion. We tested, for these runs, in 
addition to the previously used downhill simplex scheme, a new gradient-free 



optimization algorithm: the NEWUOA|49[| scheme. It is based on the con- 
struction of a higher order polynomial approximation to the function that 
needs to be optimized. 

We used a total pulse length of T = 400 a.u., larger than in the previous 
case, in order to make the nuclear movement clearly non-negligible. The sine 
series expansion contained in this case 14 components, making the parameter 
space of 13 degrees of freedom. Several initial guesses of the form fHS]) were 
tried, with uq = {3 . . . 9) -10^^ a.u.. Each initial pulse was then optimized with 
the two maximization algorithms. We observed a much faster convergence 
(roughly double) with the NEWUOA algorithm. All tests, no matter what 
maximization algorithm was used, were successful: the optimized pulse led 
to the breaking of the selected bond. We describe the results obtained for 
the case Uq = Q ■ lO^^a.u. (since all other cases showed a similar behaviour). 

The left plot of Fig. M compares both optimization algorithms. Clearly, 
the NEWUOA algorithm finds the maximum much faster. The right plot 
compares the optimized laser pulses. Here we see that the two algorithms 
found different local maxima - even if both achieved the attempted goal: the 
breaking of the selected bond. 

The left plot of Fig. [7] displays the momenta of the nuclei during the bond 
breaking test run (performed with the laser pulse of iteration step 40 of the 
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Figure 6: Left: Convergence plot for the two different optimization algorithms. DSA 
stands for downhill simplex algorithm. Right: Initial and final laser pulse, for both the 
NEWUOA and the downhill simplex algorithms. 
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Figure 7: Left: Momentum of each nucleus during the bond breaking test. Right: Nuclear 
coordinates, for the same calculation. 
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NEWUOA optimization). The right plot displays, in turn, the coordinates. 
The selected bond is broken quickly. We note that a certain degree of ioniza- 
tion occurred in all runs (the final charges oscillated between 1.1 and 2.0 a.u.) 
In all cases, most of the charge remained in the dimer fragment, permitting 
its stability. We can conclude that the inclusion of the movement of the 
nuclei in the optimization runs solves the problems found in the previous 
section, when the pulse durations are not very short. 

4-3. CG optimization for fixed nuclei 

A target constructed in terms of the momenta can also be handled with 
a gradient based algorithm, for which the QOCT equations are needed. This 
subsection describes such calculation for the same model used in the pre- 
vious two subsections. Note, however, that the QOCT equations presented 
above are valid for a quantum system, not for a mixed quantum-classical 
one. Therefore, the nuclei must be frozen during the optimization, and in 
consequence we are restricted once again to short pulses (T = 200 a.u., a case 
for which we saw that the frozen nuclei approximation is justified). 

The laser pulse was represented by the Fourier series ^2M . and further 
constraints (constant fiuence, zero average field) were then implemented as 
described in Section r2.3.1l - the parameter set is then a set of hyperspherical 
angles 6. We slightly changed the definition of the target: 

jiM = {P2m-pm) - iobi[^] -P2[^]i , (49) 

in order to have linear dependence with respect to the momenta for the two 
terms in the right hand side, since we observed that this choice usually pro- 
vides better convergence. The factor "10" can also be changed, and regulates 

28 



the weight that is placed on the minimization of the momenta difference be- 
tween those atoms that must remain bound. 

Due to the time- dependent nature of the target, Eq. [11] is now inhomo- 
geneous, and the evolution of the auxiliary wave function x is governed by 
the following equations: 

xWx.T) = 0. (51) 



The gradient VeG[9] can be calculated by Eq. (I9]).|5j] This gradient can then 
be used to perform a conjugate gradients |5o| optimization. 

We performed a number of runs with this scheme, and the results did not 
differ qualitatively of the results obtained with the gradient-free algorithm: 
partial ionization, and successful bond-breaking in about half of the runs. 
The purpose of these calculations was to make a comparison regarding the 
computational efficiency, and therefore we only show results corresponding to 
one run that was performed with identical parameters with both optimization 
schemes. 

The left plot of Fig. [8] compares the convergence for the two methods. 
The NEWUOA algorithm reached the maximum after about 60 propagations, 
whereas the CG method just needed about 25 propagations. This was typical, 
in all runs, the NEWUOA method needed about twice the computing time to 
reach convergence, (note that in the CG case, each propagation corresponds 
to either a backwards or a forwards propagation, which require roughly the 
same computer time). 

The right plot of Fig.[8]shows the optimised laser pulses for both methods. 
One can see that the two pulses look rather similar. Nevertheless, there are 
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Figure 8: Left: Comparison between the CG and NEWUOA optimization. Right: Com- 
parison of the optimized laser pulses. 

some small differences in the optimised pulses which became noticeable in 
the ionisation of the system: while the electronic charge decreased to about 
1.65 a.u. when irradiating the system with the pulse obtained with the CG 
optimisation run, we got a decrease of the electronic charge to about 1.3 a.u. 
when the NEWUOA pulse was used. 

We can conclude that a gradient-based technique such as CG is also 
applicable to this problem, and is even more efficient, despite the compli- 
cations due to the necessity of backwards propagating an inhomogeneous 
Schrodinger-like equation. Unfortunately, the scheme cannot yet be applied 
to longer pulses in which the nuclei should be allowed to move. In those 
cases, the nuclear equations of motion must then be included, as well as the 
electronic quantum equation, in the OCT formalism. Work along these lines 
is in progress. 
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Figure 9: Sketch of the 5-atomic ID test model. 

5. Selective bond breaking of ID chains 

A more stringent test on the methodology consists of attempting to obtain 
different sized fragments in longer ID atomic chains. We now show calcula- 
tions of five equal mass atom chains, for which we attempt to break the chain 
into either 4+1 or 3+2 fragments (see Fig. [9]). The chain consists of 5 Hydro- 
gen nuclei; as in previous section, they interact with the electrons through a 
soft Coulomb potential. We place four non-interacting electrons; this means 
that instead of one single wave function \E', we now have two doubly occupied 
orbitals ipi,ip2- The construction of the control target was based on the same 
ideas discussed earlier: maximising or minimising momenta differences. For 
example, for the 4+1 cleavage attempt: 

o^l[^l,^2] = (P4[^Al,^2] -P5[^l,'02]) 
3 

- 10 ^ [^1 , ^2] - Pi+1 [^1 , ^2] I , (52) 
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Figure 10: Left: Convergence plot (NEWUOA algorithm) for the 4+1 chain bond breaking 
attempt. Right: Initial laser pulse and optimized laser pulse. 

whereas for the 3+2 case: 

4 

-10 ^ \pi['ilJi,t/j2] - Pi+i[i^i,^2]\ , (53) 

i=l,i^3 

In this case, we considered moving nuclei and we apphed a gradient free 
optimization by making use of the NEWUOA algorithm. Again, we used 
fHHj) as initial pulse with a pulse duration of T = 400 a.u.. The pulse was 
represented by the constrained sine series and in this case we restricted the 
parameter search space to 11 hyperspherical angles. The following initial 
parameters have been tested: uq = {3 . . . 6) ■ 10~^ a.u.. 

For the 4+1 bond breaking attempt, almost all runs were successful (for 
two of them the field was too weak to remove any nucleus). There was no 
ionization, and all the electronic charge remained by the 4 nuclei, while one 
proton separated away. The plots in Fig. [10] and [11] correspond to the run 
with uq = 6 ■ 10~^ a.u.. The other successful runs were qualitatively similar. 
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Figure 11: Left: Momentum of each nucleus during the optimization run, for the 4+1 
bond breaking attempt. Right: Coordinates of the nuclei during the bond breaking test 
run. 

For the 3+2 bond breaking attempt, the results were different. The op- 
timization converged for all runs. However, only 2 of 10 bond breaking test 
runs were successful. In the other cases, we either obtained no ionization 
and unwanted 4+1 separation like in the previous case, or else substantial 
ionization and Coulomb explosion of the full system. 

The plots shown in Fig. [12] correspond to one successful run, namely 
that with uq = Q ■ 10~^a.u.. Fig. [13] shows the corresponding electronic 
density distribution and the coordinates of the nuclei at different times. At 
t = 300 a.u., an ionization of the system is observed. In fact, we found that 
a certain ionization was needed to remove the two nuclei. In this particular 
case, the electronic charge decreased from 4.0 a.u. to 2.3 a.u. during the laser 
pulse. This ionization necessarily implied an unwanted effect, namely that 
nucleus 4 and 5 were not bound to each other anymore after removing them 
(see plot for t = 600 a.u.). 
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Figure 12: Left: Momentum of each nucleus during one of the the 3+2 optimization runs. 
Right: Coordinates of the nuclei during the corresponding bond breaking test run. 

Therefore, for this particular choice of model, search space and algorithm, 
the optimization runs did not succeed. However, we expect that this can be 
cured in a number of ways, since there is a large freedom to be explored 
regarding the definition of the target functional. For example, the intro- 
duction of the ionization in the definition (prevention or encouragement of 
ionization) could help to avoid undesired effects caused by it. As a final re- 
mark, we mention that we performed further tests with atomic chains with 



different masses; 



54f | in those cases, it was found that the momenta should 



be substituted by the velocities in order to obtain better results. 



6. H+ 

The next example is a more realistic molecular description: a 3D calcu- 
lation for the molecule, considering interacting electrons. Fig. [H] shows 
the geometry of this molecule; jsS] it has an equilateral shape with an edge 
length of 1.64 a.u.. 
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Figure 13: Electronic density, and nuclear coordinates at different times during the bond 
breaking test run for one of the 3+2 chain cleavage attempt. The laser pulse duration was 
T = 400 a.u. 
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Figure 14: Sketch of H3". The dashed line indicates the separation plane, where the 
molecule ought to be broken. The normal vector n lies in the molecular plane and is 
perpendicular to the separation plane. The laser polarization is parallel to n, and the 
direction of the optimized momenta pi is parallel to n as well. 
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The two electrons were treated with TDDFT and the exchange-correlation 
potential was approximated by the ALDA. The motion of the nuclei was 
treated classically. The electron-nucleus interaction was described by pseudo- 
potentials - in this case, obviously, the pseudo-potentials are not used to re- 
move any core electrons, but as a means to smooth the Coulomb singularity. 

We tried to obtain a laser pulse which removes one particular nucleus, 
leaving a bound Hydrogen molecule, by making use of the same kind of 
momentum target described above. Fig [H] shows the directions in which the 
momenta were optimized; the control target is defined as: 

Jim = n-{p,m-p2m) - ip^m-p^m , (54) 

where \1/ is the Kohn-Sham orbital occupied by the two electrons. Note that 
we this functional is an explicit functional of the density. 

We used the gradient free NEWUOA algorithm for the optimization, and 
did not neglect the nuclear movement. The initial pulse was chosen to be 
in the form given by Eq. ( HHj) . and the pulse duration was T = 400 a. u.. In 
this 3D case, we also have to specify the laser polarization, which was chosen 
parallel to n. The parametrisation used to represent this laser pulses was the 
constrained sine series. 

We display in Figs. [15] and [16] the results corresponding to one typical 
optimization run, corresponding to an initial guess with ojq = 3 ■ 10^^ a. u.. 
It can be seen how the convergence is rather fast, and the obtained pulse 
cuts the molecule in the desired way. The electronic charge decreased fromn 
2.0 a.u. to around 1.5 a. u. 
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Figure 15: Left: Convergence plot for the Hg" example. Right: Initial and optimiized 
(corresponding to iteration step 40) laser pulses. 



7. CHaNH^ 



A more complex molecule is CH2NH2^, the "methaniminium cation". The 



loss of as well as from CHgNE 
both experimentally and theoretically. 



has been extensively investigated, 
SgI . Is?! Our goal was the former, the 



removal of one of the protons, the one that binds to the Nitrogen nucleus 
(this process leads to CHgNH, "methylenimine" , see Fig. [T7|l . 

We started our calculations from the ground state in which CHgNHg*" has 
a planar shape (see Fig. [T7I) . Then, the simulation of the molecule dynamics 
of CHgNHg*' was performed similarly to that of Hg^, with the described mixed 
quantum classical description on top of TDDFT. The exchange-correlation 
potential was approximated by the ALDA. The potentials of the nuclei were 
described by pseudo-potentials (in this case, this means that the two core 
electrons of C and N are frozen). 

Since we are now working with a molecule that contains nuclei with dif- 
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Figure 16: Isosurface plot of the electronic density and the corresponding positions of the 
nuclei during the bond breaking test run at different times. The isosurface was plotted at 
a density of 0.07 a.u.. The laser pulse duration was T = 400 a.u.. 
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Figure 17: Sketch of CH2NH2 . The dashed hne indicates the separation plane, where 
the molecule ought to be broken. The normal vector n lies in the molecular plane and is 
perpendicular to the separation plane. The laser polarization as well as the directions in 
which the velocities Vi were optimized are parallel to n. 

ferent masses, we will define our target in terms of the velocities, instead of 
using the momenta (the nuclear labels are defined in Fig. [T7ll : 

6 

Ji[n] = n ■ {vi[n] - V2[n]) - 10 ^ |-iT2M - Vi[n\\. (55) 

i=3 

Again, this target functional is an explicit functional of the electronic density; 
this is important conceptually since we are using TDDFT, where the many- 
body wave function is not easily accessible. In the previous equation, we 
have show explicitly this functional dependence on the density. The normal 
vector n as well as the laser polarization direction were chosen to be parallel 
to the bond axis between the Nitrogen nucleus and the Hydrogen nucleus. 

Again, we used the NEWUOA algorithm and the form given in Eq. fHHl) 
for the initial pulse; the pulse duration was T = 400 a.u. The electric field 
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Figure 18: Left: Convergence history of the CH2NH2'' dissociation attempt, for which the 
NEWUOA algorithm was used. Right Initial and optimized laser pulses for this case. 

was expanded in a sine series, and the constrained sine series parametrisation 
was used once again (this time, with 10 degrees of freedom). As usual, we 
performed optimisations with a number of initial guesses, varying frequencies 
and amplitudes (but keeping the fluence constant). Only one of the attempts 
was successful, namely that with the initial frequency = 3 ■ 10^^ a.u. The 
plots in Fig. [18] and [19] correspond to this successful run. The electronic 
charge decreased from 12.0 a.u. to 11.0 a.u. in this run. In the other cases, the 
amplitudes of the optimised electric fields were either too small or too large: 
too small electric fields merely led to oscillations of the nuclei around their 
equilibrium positions, whereas too large fields, on the other hand, caused 
high ionisation, which led to unintended dissociations. 

8. Conclusions 

This work addresses the challenge of selective photochemistry by means 
of high intensity shaped ultra-short laser pulses. The rapid experimental 
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Figure 19: Isosurface plot of the electronic density and the corresponding positions of the 
nuclei. The isosurface value of the density was 0.045 a. u.. The laser pulse duration was 
T = 400 a.u.. 
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advances in forming laser pulses of almost arbitrary shape call for reliable 
theoretical tools to predict optimal pulse shapes for certain predefined tasks. 
To achieve this goal, our basic strategy is to combine the mathematical frame- 
work of optimal control theory with a mixed quantum-classical description 
of the molecular degrees of freedom: The electronic response of the system 
is described from first principles using TDDFT while the nuclear degrees of 
freedom are governed by classical equations of motion with Ehrenfest forces 
that mediate the coulpling to the electronic degrees of freedom. 

The task that the laser pulse is supposed to perform has to be formulated 
in terms of a "control target" or "target functional" to be maximized by 
the optimal pulse. Usually a given task, like breaking a selected bond, can 
be formulated in terms of several possible target functional. This is where 
mathematical intuition and physical creativity come into play. The mixed 
quantum-classical description employed in this work lends itself to formulat- 
ing the target functional for bond breaking in terms of the classical nuclear 
degrees of freedom. We have explored target functionals based on the classi- 
cal forces acting on the nuclei, either considering their value at the end of the 
laser pulse, or considering their integrated value over the full propagation. 
This latter case means that the target functional depends on the nuclear 
momenta at the end of the pulse. The results show a clear superiority of 
the momentum-based target functional. This makes perfect sense because 
the oscillatory character of the forces makes their value at a single point in 
time less relevant than the integrated values. For molecules with different 
nuclei, it turns out to be better to define the targets in terms of the nuclear 
velocities rather than the momenta. 
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After defining the microscopic description of the system and choosing the 
control target functional, there is still ample freedom in the choice of op- 
timization algorithms. We have utilized two fundamentally different types: 
gradient-free and gradient based algorithms. The latter were found, not sur- 
prisingly, to perform better. They require, however, a more elaborate theory, 
since the gradient (or functional derivative) calculation involves the back- 
wards propagation of an auxiliary wave function which is particularly com- 
plicated when the basic equation of motion is non-linear (like in TDDFT). 33 1 

The calculations presented for 133"*" and for CH2NH^ clearly demonstrate 
that selective bond breaking can be achieved with the target functionals 
and optimization algorithms developed in this work. An immediate task 
for the future will be the application to larger molecules. Furthermore, one 
may consider the definition of refined target functionals in order to prevent 
that the removed fragments break apart later. For example, a term that 
enhances the electronic charge localization between the nuclei of the removed 
fragments could be included in the target functional. Work along these lines 
is in progress. 
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